function [rhs ind_bcD Energy_boundary] = ex33_compute_rhs(x_nodes,options,parameters,u_nodes)
% function [rhs ind_Dirichlet] = ex33_compute_rhs(x_nodes,options,parameters);
%
% This function computes the right hand side rhs
%
% Input
%    options   : lme options
%
% Output
%    rhs    : right hand side
%    ind_bcD: position of Dirichlet's nodes


% Number of nodes
nPts = length(x_nodes);
% The rhs is initialized
rhs  = zeros(2*nPts,1);
eps  = 1.0e-5;

if nargin<4
    u_nodes = zeros(nPts,2);
end


% Material parameters
a  = abs(parameters.a);
L  = parameters.L;


% The BCs nodes are classified
ids   = (1:nPts);
% (y=0)
ind_bc1 = ids(abs(x_nodes(:,2))<eps);
% (x=L)
ind_bc2 = ids(abs(x_nodes(:,1)-L)<eps);
% (y=L)
ind_bc3 = ids(abs(x_nodes(:,2)-L)<eps);
% (x=0)
ind_bc4 = ids(abs(x_nodes(:,1))<eps);


% options for 1D integration and 1D lme
opt_int1D.orderGL = 30;
opt_int1D.pruning = 0;
options1D         = options;
options1D.dim     = 1;
options1D.range_n = options.range_n;
Energy_boundary = 0;



figure(2)

%% Boundary nodes corresponding to x=L
y1D           = x_nodes(ind_bc2,2);
[y1D,ind_y1D] = sort(y1D);
ind_bc2       = ind_bc2(ind_y1D);
u1D           = u_nodes(ind_bc2,1);
v1D           = u_nodes(ind_bc2,2);

[y_s1D w_s1D] = MakeGLSamples1D(y1D, opt_int1D);
[beta_n1D range_n1D] = NodalThermalization(y1D, options1D,[1 length(y1D)]);
options1D.beta_n     = beta_n1D;
options1D.range_n    = range_n1D;
[p_lme1D dp_lme s_near1D]= ex33_wrapper_lme(y1D,y_s1D,options1D);
rhs_aux              = zeros(2*length(y1D),1);

for i=1:length(y_s1D)
    w_s              = w_s1D(i);
    i_nears          = s_near1D{i};
    p_nears          = p_lme1D{i};

    u1D_s = 0;
    v1D_s = 0;
    for k=1:length(i_nears)
        u1D_s = u1D_s + u1D(i_nears(k))* p_nears(k);
        v1D_s = v1D_s + v1D(i_nears(k))* p_nears(k); 
    end
    
    [s_xx,s_yy,s_xy] = ex33_AnalyticalStress(L,y_s1D(i),a);
    t_ = [s_xx s_xy];
    rhs_aux(2*i_nears-1) = rhs_aux(2*i_nears-1) + t_(1)*p_nears*w_s;  
    rhs_aux(2*i_nears)   = rhs_aux(2*i_nears)   + t_(2)*p_nears*w_s;  

    Energy_boundary = Energy_boundary + t_(1)*w_s*u1D_s + t_(2)*w_s*v1D_s;

end
rhs(2*ind_bc2-1) = rhs(2*ind_bc2-1) + rhs_aux(1:2:end);
rhs(2*ind_bc2)   = rhs(2*ind_bc2)   + rhs_aux(2:2:end);


% Boundary nodes corresponding to y=L
x1D           = x_nodes(ind_bc3,1);
[x1D,ind_x1D] = sort(x1D);
ind_bc3       = ind_bc3(ind_x1D);
u1D           = u_nodes(ind_bc3,1);
v1D           = u_nodes(ind_bc3,2);
[x_s1D w_s1D] = MakeGLSamples1D(x1D, opt_int1D);
[beta_n1D range_n1D] = NodalThermalization(x1D, options1D,[1 length(x1D)]);
options1D.beta_n     = beta_n1D;
options1D.range_n    = range_n1D;
[p_lme1D dp_lme s_near1D]= ex33_wrapper_lme(x1D,x_s1D,options1D);
rhs_aux              = zeros(2*length(x1D),1);
for i=1:length(x_s1D)
    w_s              = w_s1D(i);
    i_nears          = s_near1D{i};
    p_nears          = p_lme1D{i};

    u1D_s = 0;
    v1D_s = 0;
    for k=1:length(i_nears)
        u1D_s = u1D_s + u1D(i_nears(k))* p_nears(k);
        v1D_s = v1D_s + v1D(i_nears(k))* p_nears(k); 
    end
    
    [s_xx,s_yy,s_xy] = ex33_AnalyticalStress(x_s1D(i),L,a);
    t_ = [s_xy s_yy];
    rhs_aux(2*i_nears-1) = rhs_aux(2*i_nears-1) + t_(1)*p_nears*w_s;  
    rhs_aux(2*i_nears)   = rhs_aux(2*i_nears)   + t_(2)*p_nears*w_s;  

    Energy_boundary = Energy_boundary + t_(1)*w_s*u1D_s + t_(2)*w_s*v1D_s;

end
rhs(2*ind_bc3-1) = rhs(2*ind_bc3-1) + rhs_aux(1:2:end);
rhs(2*ind_bc3)   = rhs(2*ind_bc3)   + rhs_aux(2:2:end);

ind_bcD = [2*ind_bc1];
ind_bcD = [ind_bcD 2*ind_bc4-1];


